#ifndef __MATRIX_H__
#define __MATRIX_H__

#include "Vector.h"

template <int N, class Real>
class Matrix
{
public:
	// Construction. In the second constructor
	Matrix();
	Matrix(bool zero);
	Matrix(const Matrix& m);
	
	void MakeZero();
	void MakeIdentity();
	void MakeDiagonal(const Real* diag);

	// member access
	operator const Real* () const;
	operator Real* ();
	const Real* operator[] (int row) const;
	Real* operator[] (int row);
	Real operator() (int row, int col) const;
	Real& operator() (int row, const Vector<N,Real>& v);
	void SetRow(int row, const Vector<N,Real>& v);
	Vector<N,Real> GetColumn(int col) const;
	void GetColumnMajor( Real* major ) const;

	// assignment
	Matrix& operator= (const Matrix& m);

	// comparison
	bool operator== (const Matrix& m) const;
	bool operator!= (const Matrix& m) const;
	bool operator< (const Matrix& m) const;
	bool operator<= (const Matrix& m) const;
	bool operator> (const Matrix& m) const;
	bool operator>= (const Matrix& m) const;

	// arithmetic operations
	Matrix operator+ (const Matrix& m) const;
	Matrix operator- (const Matrix& m) const;
	Matrix operator* (const Matrix& m) const;
	Matrix operator* (Real scalar) const;
	Matrix operator/ (Real scalar) const;
	Matrix operator- ()const;

	// arithmetic updates
	Matrix& operator+= (const Matrix& m);
	Matrix& operator-= (const Matrix& m);
	Matrix& operator*= (Real scalar);
	Matrix& operator/= (Real scalar);

	// matrix products
	Matrix Transpose() const;
	Matrix TransposeTimes( const Matrix& m ) const;
	Matrix TimesTranspose( const Matrix& m ) const;

	// matrix-vector operations
	Vector<N,Real> operator* (const Vector<N,Real>& v) const;
	Real QForm( const Vector<N,Real>& u, const Vector<N,Real>& v) const;

protected:
	// for indexing into the 1D array of the matrix, col + n*row
	static int I ( int row, int col );

	// support for comparisons
	int CompareArrays( const Matrix& m ) const;

	Real entries[N*N];
};

template <int N, class Real>
Matrix<N,Real> operator* (Real scalar, const Matrix<N,Real>& m);

template <int N, class Real>
Vector<N,Real> operator* (const Vector<N,Real>& v,
						  const Matrix<N,Real>& m);
#endif



